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This paper is a first report on the determination of Aj^g from lattice simulations with 2+1+1 
twisted-mass dynamical flavours via the computation of the ghost-gluon coupling renormalized in the 
MOM Taylor scheme. We show this approach allows a very good control of the lattice artefacts and 
confirm the picture from previous works with quenched and N/=2 twisted-mass field configurations 
which prove the necessity to include non-perturbative power corrections in the description of the 
running. We provide with an estimate of Aym in very good agreement with experimental results. 
To our knowledge it is the first calculation with a dynamical charm quark which makes the running 
up to a s {Mz) much safer. 
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I. INTRODUCTION 



QCD, the theory for the strong interactions, can be confronted with experiments only after providing it with a few 
inputs: one mass parameter for each quark species and the only surviving parameter in the limit of massless quarks, 
namely Aqcd, the energy scale used as the typical boundary condition for the integration of the Renormalization 
Group equation for the strong coupling constant. Thus, contrary to its running that can be computed in perturbation 
theory, the value of the renormalized strong coupling at any scale, or equivalently Aqcd, has to be taken from 
experiment. 

The QCD running coupling can be also obtained from lattice computations, where the lattice spacing replaces Aqcd 
as a free parameter to be adjusted from experimental numbers: masses, decay constants, etc. Different methods have 
been used for the lattice calculation of Aqcd- Among the most extensively applied, we can enumerate the implemention 
of the Schrodiger functional scheme (see, for instance, [IH1] and references therein), those based on the perturbative 
analysis of short-distance sensitive lattice observables as the inter-quark static potential (see for instance [6]), 
heavy-quark potential, Wilson loops or small Creutz ratios expanded in the "boosted" lattice coupling (see [7HT0] and 
reference therein) or the vacuum polarization functions [111 I12j : and, in particular, those based on the study of the 
momentum behaviour of Green functions (see |13H19| and references therein). In previous studies we compared the 
behaviour of the 2-gluon and 3-gluon Green functions as a function of the renormalization scale with the perturbative 
predictions. This allowed us to get estimations for as and Aqcd; but it also revealed the presence of non-perturbative 
power corrections which we interpreted as coming from the dimension-two non-zero Landau-gauge gluon condensate 
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in an OPE approach 1 . 

In the last few years, several authors of this paper have been pursuing a program to study the running of the 
strong coupling, and so evaluate Aqcd> grounded on the lattice determination of the ghost-gluon coupling in the 
so-called MOM Taylor renormalization scheme. The main advantage of this ghost-gluon approach is that the lattice 
computation of the coupling only involves the calculation of two-point correlators, which yields a very good control 
of the lattice artefacts over a large momentum window, mainly owing to the i?(4)-extrapolation prescription [23"]l24j. 
and then for a precise checking of the running. We have first analyzed the pure Yang-Mills case (N/ = 0) [55] 
and next extended the analysis to the case in which twisted Ny = 2 dynamical quarks were included in the lattice 
simulations [2"TH2"§] . Now, for the first time, we apply the same approach (outlined in sec. [TTJ) to study the strong 
coupling by dealing with lattice simulations with two light degenerate twisted-mass flavours and a heavy doublet 
to include the strange and charm dynamical quarks. This is done within the framework of the European Twisted 
Mass (ETM) collaboration from where we used several ensembles of gauge fields for different bare lattice couplings, 
twisted masses and volume to conclude that: (i) the running description of the data requires to take into account non- 
perturbative power corrections, which appear to behave as OPE [30l [31] predict when a non- vanishing landau-gauge 
dimension-two gluon condensate is considered; and (ii) only after taking into account the non-perturbative power 



corrections, the lattice estimate for Aqcd strikingly agrees with the experimental result, as can be seen in sec. IV 



II. ABOUT THE PROCEDURE 



We shall follow the procedure described in detail in refs. [551 HI] to extract an estimate of A^g from the non- 
perturbative lattice determination of a,T(q 2 ), the running strong coupling renormalized in the MOM Taylor scheme 
and Landau gauge. Let us recall briefly how the procedure works. The Taylor coupling is defined by 

a T (^) ee g ^fl = lim ^G( M 2 , A 2 )F 2 ( M 2 , A 2 ) , (1) 

47T A->-oo 47T 

where F and G stand for the bare ghost and gluon dressing functions in Landau-gauge. As was thoroughly explained 
in ref. [55] (see also the appendix A of ref. [32 ), the well-known Taylor's paper [33] proved that, at any order in 
perturbation, the proper ghost-gluon vertex trivially takes its tree-level form when the incoming ghost momentum 
vanishes, this implying that the renormalisation constant for this proper ghost-gluon vertex is just equal to 1 in the 
MOM-like scheme defined by the particular kinematics with a zero-momentum incoming ghost, i.e. Taylor scheme. 
This is not only true in perturbation but it can be also concluded that taking the limit of a vanishing incoming ghost 
momentum drops any non-perturbative correction away from the whole proper ghost-gluon vertex, as was discussed 
in ref. |34j . Thus, Eq. (JlJ for the Landau-gauge Taylor-scheme running coupling can be straighforwardly derived from 
this last result. 

The ghost and gluon dressing functions will be here obtained from N/=2+l+l gauge configurations for several bare 
couplings, light twisted masses and volumes. Contrarily to the analysis performed in ref. |27j . the interplay of light 
and heavy quark mass and UV cut-off effects makes a chiral extrapolation harder. Further studies arc underway for 
a better control of this point. Thus, for the scope of this paper, we will content ourselves with an estimation of the 
uncertainty due to the quark mass effects. 

It should be emphasized that a crucial role is played by the appropriate elimination of discretization artefacts to 
provide us with reliable and exploitable results. A first step consists in curing the artefacts which are due to the 
breaking of the rotational invariance on the lattice, where the remaining symmetry is restricted to the H(4) isometry 
group. For this purpose, we perform the so-called i/(4)-extrapolation procedure [53J[5U[35] that leaves us with 



/ n M \ f) n hatt 



a 2 ^ + ■•■ (2) 

DM „ 



where p' 4 ' = YliPi 1S the first iJ(4)-invariant (and the only one indeed relevant in our analysis). Thus, we first average 
over any combination of momenta being invariant under i?(4) (H(4) orbit) and extrapolate then to the "continuum 



The possible phenomenological implications in the gauge-invariant world of such a dimension-two Landau-gauge gluon condensate and 
in connection with confinement scenarios has been also largely investigated, as can be seen for instance in ref. |20l |21| . This condensate 
has been also related to the QCD vacuum properties through the instantons liquid model |22| 
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case" , where the effect of a 2 p^ must vanish, by applying Eq. ^ for all the orbits sharing the same value of p 2 , with the 
only assumption that the slope depends smoothly on a 2 p 2 and can be fitted to a polynomial form from the whole set of 
lattice data. Furthermore, the H (4)-artefact-free lattice coupling, o+(a 2 p 2 ) might differ from the continuum coupling 
by some 0(4)-invariant artefacts, as shown for example in the lattice analysis of the quark-propagator renormalization 
constant [3"BH3"5] . This leads us finally to write: 

a T (a 2 p 2 ) = a T {p 2 ) + c a2p2 a 2 p 2 + <D{a A ) , (3) 

where c a 2 P 2 should be fitted from the lattice data and verify the appropriate scaling from the simulations with different 
bare couplings, ft; while aj- is the lattice prediction to be compared with the continuum OPE formula for the Taylor 
strong coupling [2"7] , 



a T (p 2 ) 




1 7 " /fh gfrtiKA*)^ 



4(/V, 



2 



(4) 



where 7^ can be taken from [39l |40] to give for N/ = 4, 



A 2 in _ 27 27 

1-70 132-8^ - 100 ; (5) 

and, applying the same method outlined in the appendix of ref. [27] . one can take advantage of the O(a 4 )-computations 
for the Wilson coefficients in ref. [40] . and obtains 



R (a, a ) = (1 + 1.18692a + 1.45026a 2 + 2.44980a 3 ) (l - 0.54994a - 0.13349a 2 - 0.10955a 3 ,) , (6) 
for q = 10 GeV. The purely perturbative running in Eq. @ is given up to four-loops by [JT] 



pot \ p 2 t p*t 2 U &w 2 J Pi 4 



pert, 2 s ^ I, Pi log(t) , PI 1 {(, M lY , P2P0 5 

* T ^''Poty PI t + J 
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' ' A ' UPA ' 21og 3 (t) + 51og 2 (t)+ 4-6^° log(t)-l 



( 7 ) 



(Pot)* 2ft 2 V/V V \ Pi 



with t = In jt and the coefficients of the p- function in Taylor-scheme [251 H2] • As for the Aqcd parameters in 
Taylor-scheme and MS, they are related through [27] 

507 - 407V/ 

^"792-48/V/ = . 5608 3 2 ; (8) 

for the N f = 4 case 2 . Thus, only three parameters, g 2 (A 2 ), A^g and the coefficient for the <9(4)-invariant artefacts 
c a 2p2, remain free to be fitted through the comparison of the prediction given by Eqs. ( 3|4 | and the lattice estimate 
of Taylor coupling after 7J(4)-extrapolation. 



III. THE LATTICE SET-UP 



As already mentioned, we obtain a^ att by Eq. jl]) from the ghost and gluon propagators computed from the 
gauge configurations simulated at several lattices with Ny =2+1+1 mass-twisted lattice flavors [13] by the ETM 



2 It should be noted that, although Ay is gauge-dependent, A^jg is not. Of course, the conversion factor of both parameters to each other 
is also gauge-dependent. Then, as far as this conversion factor can be exactly determined because of the RGE invariance of Aqcd, the 
choice of any gauge for the lattice determination of Ay (Landau gauge in our case) is irrespective for the final determination of Aj^r . 
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collaboration jlUHS]. I n the gauge sector, we use the Iwasaki action and compute the propagators as described in 
refs. |27j . while for the fermion action we have 

Si = a 4 ^2xi( x ) (d w [U] + m 0t i + ifxij 5 r^j xi(%) ( 9 ) 

for the doublet of degenerate light quarks [46] and 

Sh = a A ^Xh{ x ) [d w [U] + mo M + i/i^Ti + Hsr^j Xh{x) (10) 

for the heavy doublet. D\y[U] is the standard massless Wilson Dirac operator. The lattice parameters for the 
ensembles of gauge configurations we used are given in tab. [TJ Tuning to maximal twist is achieved by choosing a 
parity odd operator and determine n cr i t such that this operator has a vanishing expectation value. One appropriate 
quantity is the PCAC light quark mass and we demand mpcAC = 0. We refer the interested reader to refs. [44l [45] 
for more details about the set-up of the twisted mass lattice simulations. 





Werit 






a/As 


{L/af x T/a 


confs. 


1.95 


0.1612400 


0.0035 


0.135 


0.170 


32 3 x 64 


50 




0.1612400 


0.0035 






48 3 x 96 


40 




0.1612360 


0.0055 






32 3 x 64 


50 


2.1 


0.1563570 


0.0020 


0.120 


0.1385 


48 3 x 96 


40 



TABLE I: Lattice set-up parameters for the ensembles we used in this paper. They correspond with the ones coded as B35.48, 
B35.32, B55.32 and D20.48 in Tab. 1 of ref. [45] . The last column stands for the number of gauge field configurations we used. 



IV. THE RESULTS OF THE ANALYSIS 



A. Curing the f/(4)-artefacts 



The first stage of the analysis, as explained above, consists in the application of the if(4)-extrapolation to cure 
the main type of discretization artefacts, namely the ones coming from the breaking of the rotational symmetry. 
These effects appear to be very visible in Fig. [l]a, where we plot the ghost dressing function before and after H(4)- 
cxtrapolation in terms of the square of the momentum in lattice units. The classical "fishbone" structures generated by 
the different H (4) orbits corresponding to the same continuum momentum can be strikingly seen before applying the 
extrapolation. Moreover, it is obvious that, had we rather applied some sort of average over a "democratic" selection 
of the orbits at all physical momenta, the resulting ghost dressing function would have shown an anomalously flat 
behaviour, with no indication of the perturbative logarithm, in the large momentum region. All these anomalies 
appear to be strikingly cured by the i/(4)-extrapolation prescription. It should be noted that a very important input 
to apply properly such a prescription comes from the many orbits at our disposal on a large momenta window (this 
implies to Fourier transform over large momenta). It is beneficial because, for the same price, the output is a very 
clean signal over a large momenta window that permits a very precise checking of the running behaviour. The results 
of the extrapolation for the coupling in Taylor scheme are plotted against the momentum in lattice units in Fig. [I] (b) 
for the four lattice ensembles analysed. 



B. Flavour mass, finite volume and 0(4)-invariant artefacts 



After ff (4)-extrapolation, a main part of the discretization artefacts are supposed to be removed and, if so, all the 
curves for the running coupling in Fig. [l](b) should appear superimposed after being re-expressed in terms of the 
momentum in physical units, at least over the momentum region where other artefacts do not play any significant 
role. In particular, with no further conversion, the three ones for /3 = 1.95 should coincide with each other (and they 
are close to) and the one for f3 = 2.1 should also after the appropriate rescaling of the momentum in coordinates 
axis. On the other hand, the comparison of the running for the three ensembles with same j3 parameter shows that 
no relevant finite volume effect happens above a lattice momentum of the order of 0.5, while a small flavour mass 
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F(p 2 ) qt(p 2 ) 




6 8.5 1 1,5 2 2*5 3 3.5 A 4.5 G 0.5 1 1.5 2 2.5 



a\p)p 2 a(P)p 
(a) (b) 

FIG. 1: (a) An example of the ghost dressing function lattice data (case: j3 = 1.95, afii = 0.035, L = 48) before (black points) 
and after (red points) H (4)-extrapolation, plotted in terms of the square of momentum in lattice units, (b) The running of the 
coupling as a function of the momentum in lattice units for the four ensembles of lattice data: /3 = 1.95 with a/i; = 0.0035 at 
48 3 x 96 (red) and 32 3 x 64 lattices (blue), with am = 0.0055 at 32 3 x 64 (violet) and f3 = 2.1 with a/i; = 0.0020 at 48 3 x 96 
(green). 



effect appears to be visible: the two ensembles with same bare mass (a/i/ = 0.0035) and different volumes appear 
nicely superimposed while the ones with same lattice volume (32 3 x 64) and different bare masses seem to require an 
additional fine rescaling. Then, we apply a rescaling factor to the lattice momentum, when needed, to render the four 
curves coincident and show the results in Fig. [2j We choose to rescale all the data to those of the ensembles with 
a/ii = 0.0035 at j3 — 1.95, since their lattice spacing seems rather safely established and we have two different volumes 
which agree fairly well. We thus obtain an optimal rescaling factor of 1.07 for lattice momenta with a/i/ = 0.0055 at 
P = 1.95 and a factor of 1.36 (containing both the ratio of lattice spacings and possible flavour bare mass effects) for 
those with am — 0.0020 at j3 = 2.1. The agreement is indeed impressive and make us conclude that the flavour bare 
mass effect can be absorbed by physical calibration of the lattice spacing and be either removed by chiral extrapolation, 
when possible, or included in the calibration systematic uncertainty. Nevertheless, it is also manifestly shown by large 
lattice momentum pattern of data (see the right plot of Fig. |2| that the 0(4)-invariant artefacts, which cannot be of 
course cured by the if(4)-extrapolation, still survive and demand some treatment for a precise analysis of the running. 
We will proceed to remove the remaining discretization artefacts for all our lattice data sets by applying Eq. (|3| with 
the requirement that the coefficient c a 2 P 2, the correction being a lattice artefact, should be universal. This will be 
explained in the next subsection. 



C. The physical running, the gluon condensate and A 



MS 



As explained in Sec. |TTJ Eqs. (3j4| can be directly applied to fit the lattice data plotted in Fig. [T](b) with only 



three free parameters. In lattice units they are on one hand A-^a((3) and g 2 (A 2 )a 2 (/3), which depend on the lattice 
spacing at each simulation and on the physical values for Aj^g and for the Landau-gauge gluon condensate, and on 
the other hand c a 2p2 which should be the same number for any lattice data set. In the following, we will only analyze 
the two ensembles at j3 — 1.95 simulated for 32 3 x 64 lattices and the one at /3 = 2.1 for a 48 4 x 96 lattice, all of them 
sharing approximately the same lattice volume in physical units. The ensemble at f3 = 1.95 for a 48 4 x 96 lattice, 
for which we exploited a smaller number of gauge configurations and have larger statistical errors, has been used to 
check finite-size effects. Then, as we showed that no visible finite-size effect survives above a(fl)p ~ 0.5 and that the 
flavour bare mass effects can be fairly well described by a lattice calibration, we fit independently any ensemble of 
data at both j3 — 2.1 and (3 — 1.95, by imposing the coefficient c a 2 P 2 to be universal, and including a (fitted) rescaling 
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FIG. 2: The same shown in Fig.[T](b) but after the appropriate rescaling of the x-axis to superimpose the four curves as much 
as possible. The right plot is a zoom for the large momentum region to exhibit the discrepancies indicating that 0(4)-invariant 
artefacts are present. 



factor to bring the heavier mass data to superimpose with the lighter ones at /3 = 1.95. Thus, we obtain the results 
of Tab. [Tl] for the best-fit parameters. Using these values with the appropriate rescaling, we plot in Fig. [3] the fitted 
running coupling obtained for the different lattices after removing all the discretization artefacts. 







a(£) A MR 


a(P) 2 g 2 {A 2 ) 


aQ3)/a(2.1) 


C a 2p2 


1.95 


0.0035 
0.0055 


0.126(12) 
0.117(8) 


0.7(3) 
0.57(22) 


1.36(16) 
1.27(10) 


-0.0047(12) 


2.1 


0.0020 


0.092(5) 


0.40(9) 


1 



TABLE II: Best-fit parameters for the confrontation of the different lattice data and Eqs. ( 3j4 1 , as explained in the text. The 
local operator A 2 has been renormalized at qoa(2.1) ~ 3 which corresponds to go — 10 GeV when the appropriate conversion 
is applied. The errors have been computed by using a jackknife procedure. 



The universality of the coefficient c a 2 P 2 and the nature of the remaining discretization artefacts that could be seen 
in Fig. [2] can be directly checked on the data. Actually we have 



"Lt' 1 (a(l-95)p) 



4: t \ 95 (a(1.95)p) 



a 2 (2.1) 
a 2 (1.95) 



1) c a2p2 a 2 (1.95)p 2 + o(a 2 (1.95)p 2 ) 



(11) 



which means that, after an appropriate rescaling implying for all the momenta to be written in units of the lattice 
spacing at a given reference (3 (f3 — 1.95 and a/i; = 0.0035 has been choosen here), the difference between the data 
for the coupling obtained at /3 = 2.1 and those at /? = 1.95 has to be proportional to c a 2 P 2 times the square of the 
momentum when the momentum is not too large. This is illustrated in the upper-left plot in FigJSJ where the linear 
behaviour on a 2 p 2 is verified (for the two lighter masses). The slope given by the r.h.s. of Eq. (11) with the best-fit 
parameters of Tab. [TT| perfectly agrees with the data. In order to perform the subtraction in the l.h.s of Eq. ( 11 ) some 
interpolation procedure is required to estimate the values for the momenta at f3 — 2.1 from the data at (3 = 1.95. 
To this purpose, we used the expression given in Eqs. (3]4) with the best fit parameters to analytically represent the 
data at (3 — 1.95. The result for the subtraction is shown by the black points in the plot. This procedure does not 
work below a(1.95)p ~ 1.5, as the expressions Eqs. (3j4| do not represent properly the data in this region. Any other 
interpolating formula, as far as it fits well, would provide us with equivalent results ; we illustrate this point in the 
plot by using also a polynomial of fourth degree to describe the lattice data at f3 — 1.95. The result for the subtraction 
in this case is given by the blue points. The only quantity we need to estimate the coefficient c a2p 2 from the data is 
the ratio of lattice spacings. The large uncertainty given for these quantities in Tab. In] is partially a consequence of 
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our present analysis, which uses only the momentum window where the OPE prediction given by Eq. Q appears to 
be in order. As far as any UV cut-off contribution should vanish at the infinite cut-off limit for the running coupling 
defined by Eq. ([!]), all the data, properly corrected for lattice artefacts, from different lattice simulations should scale 
when expressed in terms of physical units. Had we only then be interested in obtaining the ratio of lattice spacings, 
the matching would be performed over a much larger momenta window 3 and the uncertainties would be drastically 
reduced. Anyhow, in the analysis of the present paper, we will follow a different fitting strategy, as will be seen in 
the next subsection. 





FIG. 3: (Upper-left) Check of Eq. ^TTJ from the lattice data at /3 = 1.95 {am = 0.0035) and P = 2.1 (afii = 0.0020), as explained 
in the text. (Upper-right) deviation from the lattice data free of discretization artefacts with respect to the prediction of the 
four- loop perturbative theory, with a A^jg taken from ( 14 1, plotted in terms of the perturbative running; the solid line shows the 
leading non-perturbative OPE prediction, Eq. (|15[)'s r.h.s. (Bottom-left) The departure of lattice data from the leading non- 
perturbative OPE prediction for the running coupling plotted in logarithmic scales, in terms of the momentum in lattice units of 
a _1 (1.95): a next-to-leading 1/p 6 behaviour is strikingly manifest. (Bottom-right) The physical running of the strong coupling 



obtained from the lattice data free of discretization artefacts, expressed in terms of the momentum in units of a 1 (1.95); 
solid line stands here for the best-fit with Eqs. |3|4l, while the dotted one is for the four-loop perturbative prediction. 



the 



One can perfom the matching all over the region where finite-size effects appear not to be visible, although the matching procedure 
would imply some sort of practical fit of the lattice data inside such a region. This will be the object of a forthcoming work. 



D. The global fit 



As wc concluded in the previous subsection, the analysis of the three different lattice data sets clearly indicates that 
the flavour bare mass effects can be fairly well described by a lattice calibration. This means that we can suppose 
that the lattice spacing for any bare coupling and flavour mass can be written as 

a{0,m) = a(0,O) (1+ c ow a(/3,0) 2 /z? +o(a 2 tf)) , (12) 

where c afi gives the slope for the chiral behaviour of the lattice spacing. Of course, this light-quark bare mass 
dependence for the lattice spacing must be transferred to any physical quantity like the Taylor coupling, after its 
lattice artefacts have been removed, 



OLT 



(a 2 (/3, W > 2 ) = «t(« 2 (/3,0)p 2 ) + 2a 2 (/3,0)p 2 e Q/i 



a 2 (/3,0) M 2 + 

x=a 2 (0,O)p 2 



dx 

(a 2 (/3,0)p 2 ) + i? (a 2 (/3,0)p 2 ) a 2 (/3,0V 2 , (13) 



where we also assumed the strong coupling not to "feel" any additional light-quark bare mass effect, as it is clearly 
suggested by the results of the previous subsection. In the analysis of lattice configurations for N f — 2 twisted-mass 



flavours in ref. |27j . Eq. (13) was successfully a ppli ed to extrapolate down to zero light-quark mass all the data within 
a narrow momentum window where Rq in Eq. ( |13[ ) was shown to be well approximated by a constant 4 . Here, we will 
procceed otherwise: we will take the ratios of lattice spacings, a(/3, a/ij)/a(1.95, 0.0035), from the previous subsection 
analysis and express the Taylor coupling from the three ensembles of lattice data in terms of the momentum in units 
of a(1.95, 0.0035) -1 . Then, we make a global fit for the three ensembles (for all momenta above a(1.95, 0.0035) p=1.5) 
and obtain the following results: 



a(1.95,0.0035) = 0.125(5), a 2 (1.95, 0.0035) g 2 (q 2 )(A 2 ) R . q 2 = 0.70(6), c a2p2 = -0.0046(7), (14) 

whith a best \ 2 = 103.8 for 317 degrees of freedom. Some plots resulting from the global fit can be seen in Fig. [3] 

The upper-right plot of Fig. [3] shows the lattice data, after the subtraction of the perturbative running and the 
0(4)-invariant artefacts, multiplied then by the square of the momentum in units of a -1 (1.95) and plotted in terms 
of the four-loop perturbative value of the coupling at the same momentum. To obtain the perturbative coupling, 
we apply the value of Aj^g obtained from ( 14 ) which, after being converted to physical units with lattice spacing 
taken from ref. |44] j35| , appears to be in very good agreement with the experimental result (see below) . According 
to Eq. Q, one would obtain: 

2 / I 2\ pert/ 2\\ ^9t(%)(A )ll,ql , pert , 2 \ pert/ 2\\ pert/ 2\ f a T (p 2 )\ r\ 

P (aL T {p)-a y T (p )) = A{Nl-l) ^ T ^ ^ T ^ T ™' [ a pcrt (g§) J ' 

The solid line in the upper-right plot of Fig. [3] corresponds to the r.h.s of Eq. (151 with the value for the Landau-gauge 
gluon condensate taken from Eq. (14 1. One should notice that the departure from zero for the lattice data in the plot 
can be only explained by non-perturbative contributions. Furthermore, the Wilson coefficient for the Landau-gauge 
gluon condensate in the OPE expansion successfully accounts for the non-flat behaviour from the lattice data in the 
small coupling regime. This provides with a striking indication that (and where) the OPE analysis is in order. Next, in 
the bottom- left plot of Fig.[3j we show the departure of the lattice data from the prediction given by Eqs. (3j4), plotted 



in terms of the momentum in units of a _1 (1.95), with logarithmic scales for both axes. The data seem to indicate that 
the next-to-leading non-perturbative correction is highly dominated by an 1/p 6 term. This is just a factual statement 
which might suggest either that the 1/p 4 OPE contributions are negligible when compared with the 1/p 6 ones or that 
the product of the leading 1/p 4 terms and the involved Wilson coefficients leave with an effective 1/p 6 behaviour. 
One might also guess that a different non-perturbative mechanism dominates over the momentum window where the 
dimension-four OPE condensates had to be visible. Finally, the bottom-right plot shows the physical running of the 



coupling and how well the OPE formula with the results from Eq. ( 14 ) fits the data 



4 This was found to happen for 1 < ap < 1.5, where the decreasing of the derivative of ax compensated the increasing due to the factor 
a 2 p 2 . It should be furthermore noticed that the derivative is negative, while c a(J appears to be positive from the chiral extrapolation of 
the Sommer parameter, rg/a, in ref. 1441 . This agrees with the sign of the chiral slope for the Taylor coupling in ref. 1271 . 
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E. Conversion to physical units 



The two main purposes of this paper are to show the impact of the OPE power corrections in describing the running 
of the strong coupling and to give an estimate for A^jg from lattice QCD simulations with a dynamical charm quark. 
The results for both goals can be summarized by the conversion to physical units for A^jg and the Landau-gauge 
gluon condensate. To this purpose, we will apply 

0(1.95,0.0035) = a(l f 1 5 Q "-° n °f } a(1.95,0) (16) 
a(1.95, 0) 

where the absolute calibration of the lattice spacing at /3 = 1.95, after the chiral extrapolation for the light quark 
mass, will be taken from refs. [HH2S]: a(1.95, 0) = 0.0779(2) fm; and where we approximate 

O (1.95,0.0035) _ +0 .o3 

0(1.95,0) " -° ' [ ' 

The systematic error quoted here has been estimated from the chiral extrapolation of the Sommer parameter in 
ref. [H]. There, as can be seen in plot 6.(b), one gets 

(rn/n) afJ - l=a 

, , \ nnrm ~ 1 = 0.015(11) . (18) 
(ro/a) ^ 1 ^ - 0035 

Then, if the string tension for the static interquark potential is supposed not to depend very much on the light quark 



mass, Eq. (18) gives the conservative systematic uncertainty for the deviation from 1 in Eq. (17). Thus, we apply 
Eqs. pSfll8| ) into Eq. ^ and obtain: 



A^4~ 4 = 316 ± 13 ± 8±g MeV 

Mb y 

9 2 {ll){A 2 ) R , q i = 4.5 ±0.4 ±0.2318.3 GeV 2 ; (19) 

where the first quoted error is statistical, the second one reflects the present uncertainty on the absolute calibration 
of the lattice spacing that we roughly (and conservatively) estimate to be of ±2.5 %; and the third one is for the 
chiral extrapolation of the light quark mass. More precise estimates for these systematics uncertainties will be 
accessible with more data (more simulations at different j3's and for more light quark masses). Finally, the value 
for A^ -4 and the four-loop perturbative running with the appropriate crossing of the bottom mass threshold at 
mb(mb) — 4.19lg 8 GeV [41 can be used to estimate the value of the coupling at the Z° mass, 

a s (M z o) = 0.1198(9) (5)1° ; (20) 

where the errors have been properly propagated. This is a first result that will be refined, mainly by improving the 
precision for the estimates of systematic uncertainties. However, it appears to be pretty compatible with the last 
world average given by PDG [31]: 0.1184(7). Although a more detailed comparison of our result with PDG average 
and a discussion of its implications will be left for a phenomenologically targeted forthcoming letter, we should remark 
that our result including strange and charm dynamical quarks (Np =2+1+1) appears to be slightly larger than the 
lattice estimate for Ny=2+1 staggered fermions, applied to obtain the PDG average: a>s (M z o) = 0.1183(8) |10) . 
Assuming no systematic effect to appear from the different fermion actions, the meaning for the 1-cr discrepancy of 
both central values, if any, can be explained from the procedure applied to cross the threshold from N/=2+l to 
N/=2+l+l flavours. That procedure is very well established and controlled in perturbation theory [4TJ [47] , but some 
non-perturbative effects may still appear at the charm quark running mass. 

Indeed, if we compare the result of A^ for N/ = 2 in ref. [57], the central value ranging from 310 to 330 MeV 
(depending on the lattice size calibration at different values of j3), with that for N/=2+l+l in this paper, it can 
be concluded that the effect for the running coupling of crossing the strange or charm quark thresholds is not very 
significant 5 on A^jg. On the other hand, applying just the perturbative recipe to cross the charm quark threshold 
would result in a stronger decreasing of Aj^g from N/=2+l to N/=2+l+l flavours. Thus, an enhancement for the 
estimate of as (Mzo) obtained from the N/=2+l+l lattice result for A^jg- can be understood when comparing to the 
one obtained by applying that perturbative recipe to cross the charm quark threshold with the N/=2+l lattice result. 



5 In Sec. 4.5 of ref. |27| , a significant difference between quenched (Nj = 0) and = 2 lattice results for A^g- (once lattice spacing 
is calibrated from /„■) was observed. In that case, simulations for infinite (quenching) and vanishing (chiral limit) quark flavours were 
compared, although it is well known that the chiral lmit for the quenched case is wrong. On the contrary, when comparing N/=2+l+l 
and N^=2 (or N^=2+l) cases, we deal with infinitely massive flavours (s and c) in the latter and heavy (c) or mid-heavy (s) in the 
former. Then, not to see the same significant effect cannot be too surprising. 
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V. CONCLUSIONS 

We used lattice gauge field configurations with two degenerate light and one heavy doublet of twisted mass flavors, 
produced within the framework of ETM collaboration, to compute the running strong coupling in the MOM Taylor 
scheme. In this particular renormalization scheme, the lattice computation of this coupling has the very nice feature 
of involving only propagators. This allows for a very precise control of the lattice artefacts and other systematic 
uncertainties. In particular, the so-called H (4)-extrapolation procedure which has already been proved to be very 
effective at eliminating the discretization artefacts of two-point correlation functions is efficiently at work in this 
analysis. The dominant 0(4) artefact is also rather easily isolated and eliminated. On the other hand, the renormalized 
coupling, defined in MOM taylor scheme by the combination of the ghost and gluon bare propagators, must only 
depend on the UV cut-off through residual contributions vanishing at the infinite cut-off limit. Thus, the Taylor 
coupling computed from the lattice must join the continuum prediction at this infinite cut-off limit. This is also 
a strong point to cure properly lattice artefacts and get reliable results. Thus, we evaluated the running strong 
coupling over a rather large window of lattice momenta and, after the appropriate relative calibration, confront it 
with the perturbative prediction available at the four-loop level. We clearly demonstrate the necessity to include 
non-perturbative power corrections to get an accurate description for the behaviour of the coupling constant. We 
then show these corrections follow the OPE predictions when a non-vanishing dimension-two Landau-gauge gluon 
condensate is present. This is to our knowledge the first time a Wilson coefficient is directly confronted to numerical 
results. Higher order terms are visible and seem to be dominated by 1/p 6 contributions instead of the expected 1/p 4 
ones. 

The precise comparison of lattice estimates with continuum formula allows the estimate of both the gluon condensate 
and Ajyjg. The latter is in a very good agreement with the world average of its experimental determinations provided 
by PDG, confirming that the picture we advanced, first by the analysis of quenched lattice data and next by studying 
lattice simulations with two twisted-mass flavours, agrees very well with the real world when a full QCD analysis is 
performed including a heavy doublet for strange and charm quarks. Going the other way around, with the experimental 
value for Aj^g- as an input, this approach could be used to provide quite good absolute determinations of the lattice 
spacings. Relative measurements (ratio of lattice spacings) through the superimposition of the different curves for 
a(a 2 p 2 ) over a large window can be also obtained with a good level of precision. The perturbative running from 
our results up to the Z boson mass is well known in perturbation theory and we apply the standard formula |41) to 
cross the b quark threshold. We do not need to consider the charm quark threshold since, for the first time to our 
knowledge, we have used dynamical charm in our computation. This is a very significant gain : the crossing of the 
charm threshold using perturbative QCD is quite questionable, since, as we have seen, non-perturbative effects ( OPE 
power terms) are sizeable at this energy. 

Further works are in progress. In particular data with different masses at Nf =2+1+1 and from simulations with 
Nj = 4 light flavours will soon be available to help to improve on the question of the mass effect. 
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